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Abstract 

We formulate a finite-difference time-domain (FDTD) approach to simulate electromagnetic 
wave scattering from scatterers embedded in layered dielectric or dispersive media. At the heart 
of our approach is a derivation of an equivalent one-dimensional wave propagation equation for 
dispersive media characterized by a linear sum of Debye-, Drude- and Lorentz-type poles. The 
derivation is followed by a detailed discussion of the simulation setup and numerical issues. The 
developed methodology is tested by comparison with analytical reflection and transmission coef- 
ficients for scattering from a slab, illustrating good convergence behavior. The case of scattering 
from a sub-wavelength slit in a dispersive thin film is explored to demonstrate the applicability of 
our formulation to time- and incident angle-dependent analysis of surface waves generated by an 
obliquely incident plane wave. 
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I. INTRODUCTION 



The study of obliquely incident plane wave upon planar interfaces is of fundemantal 
interest to electromagnetic (EM) wave propagation. It underlies Snell's law of refraction and 
leads to important concepts such as total reflection and Brewster's angle.- 1 ^ One can easily 
relate to the concept of an obliquely incident plane wave by the daily experience of looking 
into a mirror. In practice, oblique incidence is vastly applied in EM related applications 
such as fiber optics,- underground object detection,- and RF-human body interaction.- 

In the growing field of nanoplasmonics,-^ oblique incidence finds applications particularly 
in exciting surface plasmon polaritons (SPPs), exemplified by the common experimental 
setup in which subwavelength defects or attenuated total reflection are utilized to couple 
the obliquely incident plane wave into propagating SPPs.- 1 ^ By taking advantage of the 
incident angle degree of freedom, several experiments have demonstrated SPP near-field 
manipulation,— ~— which has been proposed as a direct approach to measuring SPP gener- 
ation efficiency— Most recently, it has been shown that SPP's can be directly generated on 
a planar metal surface by interfering incoming light beams with different incident angles in 
a four- wave mixing scheme.— 

The effects of light incident at oblique angle on sub- wavelength defects in metallic layered 
media have been studied by frequency-domain calculations based on either coupled wave 
analysis or a semi-analytical model. These references have explored the obliquely incident 
light transmission through a single defect,— ~— and the SPP generation efficiency.— 1 ^ Our 
work is largely motivated by recent experimental study of SPP dynamics excited or controlled 
by a femto-second (fs) laser pulse obliquely incident on a SPP propagating interface.— 1 ^ To 
describe such experiments, a time-domain method is desirable because of the ultrafast nature 
of the exciting or controlling laser pulse. The major challenge in developing such a method, 
is to accurately treat the oblique incidence as well as the material dispersiveness. This poses 
a special challenge for mesh-based propagation method (such as the finite-difference time 
domain method) because even in the absence of the inhomogeneous media the wave front 
not aligned with the Cartesian mesh is required to be uniform and to have arbitrary incident 
angle and time profile. 

In this paper, we develop a numerical method to rigorously treat obliquely incident plane 
wave scattering at embedded scatterers in layered dielectric and dispersive media. To the 
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best of our knowledge, such a method was not published as yet.— Targetted mainly at 
time-domain studies of EM wave phenomena that involve SPP excitation and propagation 
in metallic films, the developed method is formulated within the framework of the finite- 
difference time-domain (FDTD) method. This method has enjoyed a wide range of appli- 
cations in the field of nanoplasmonics,— ^ and its time-domain nature makes it particularly 
well suited to ultrafast phenomena. Our treatment of the oblique plane wave is an extension 
of the total field / scatter field (TF/SF) technique to describe media characterized by the 
combination of Debye-, Drude- and Lorentz-type poles.— >2£ The TF/SF technique has been 
applied successfully in the FDTD study of free-standing scatterers, layered dielectrics and 
dispersive media describable by a single Debye pole.— >2£~— It is based on the linearity of 
Maxwell's equations and decomposes the total field into an incident and a scattered field 
components,— 

By setting up an artificial boundary between the TF and SF regions in the FDTD simula- 
tions, a plane wave of arbitrary time profile and incident angle can be achieved by matching 
the known incident field at the TF and SF boundary. In presenting our method, we will 
focus on the derivation of an equivalent one- dimensional (ID) wave equations for the TF/SF 
boundary condition, suitable for various types of material dispersiveness, and explain in de- 
tail the numerical considerations involved. This will be followed by extensive numerical tests 
of the convergence properties of the method. For clarity, several important concepts from 
the previous literature are reemphasized. 

This paper is organized as follows: In Section [TTJ we derive the equivalent ID wave 
equations, show the numerical flow chart for matching the TF/SF boundary condition, 
and discuss several practical simulation details, including stability, interface treatment, and 
leakage. Section II I II tests our approach by comparison of numerical with analytical results 
for model problems. Finally, concluding remarks are provided in Section IIVI 

II. THEORY AND NUMERICAL METHOD 

In the following, we provide the equations and numerical method for solving the trans- 
verse magnetic (TM) mode in two dimensions (magnetic field perpendicular to the two- 
dimensional plane). Special emphasis is placed on the TM mode because of its relevance 
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to SPP excitation.— The numerical approach for solving the transverse electric (TE) mode 
equations is similar to that for the TM mode, and the corresponding equations are given in 
Appendix |A] The media considered are vacuum, linear dielectric media (characterized by 
a dielectric constant), and linear dispersive media (characterized by a finite sum of Debye, 
Lorentz, and Drude types of poles). 



A. TM mode wave propagation: two-dimensional and equivalent one-dimensional 
equations 

Our starting point is Maxwell's equations in the frequency domain for the TM mode, 

dE y dE x . 
ox ay 

^1 = -iue e(uj)E x , (3) 

oy 

dH z . . 

-g— = iue e(u)E y , (4) 

where the coordinate system is defined in Fig. (TJ €q is the free space permitivity, no is the 
free space permeability, and e{uS) is the dielectric function for a dispersive media, which 
reduces to a constant for vacuum and dielectric media. 

In the case studies below, we assume a dispersive medium with a single (non-zero) Drude 
pole e m = e(oo) — uj 2 d /{uj 2 + iT D u) and provide two separate sets of equations for solving 
Eqs. ([2111]). The first set of equations is based on the auxiliary differential equation (ADE) 
approach with polarization currents to account for the dispersiveness. In this case, we 
further assume that media other than vacuum are not extended into the absorbing boundary, 
which allows us to use Berenger's PML absorbing boundary condition. 31 The second set 
of equations is formulated within the general context of the Uniaxial Perfectly Matched 
Layers (UPML) absorbing boundary conditions,— and involves a different approach to treat 
the dispersiveness. In this case, we can effectively absorb the outgoing waves exiting the 
simulation domain in the dielectric and dispersive media. Separate tests have been done 
to ensure that the two approaches provide the same solution.— In the following, we assume 
that the 2D electric and magnetic fields propagate on the Yee mesh with the dependence 
«|"- = u(iAx, jAy, nAt). Ax = Ay is the size of a unit cell, and At is unit time step. 
Details on the FDTD equations in both ADE and UPML approaches are given in Appendix 



If we now consider TM mode wave propagation with obliquely incident plane wave on 
layered media with translational invariance, the 2D equations of motion can be reduced to an 
equivalent ID wave propagation problem along the direction perpendicular to the interfaces 
between the media.— >^ We proceed to derive the equivalent ID wave equation for the TM 
mode, which will serve as a means of introducing incident fields along the TF/SF boundary. 
The corresponding derivation for the TE mode is provided in Appendix [A] 

Substituting Eq. @j into Eq. (T5]) yields, 

dE * ■ tt 1 d 2 H z 

— — = -iujhqH z + . (5) 

oy iuje e{uj) ax A 

Because of the translational invariance and phase matching across the interfaces between 
different layers, d 2 H z /dx 2 = —k 2 H z , with k x being a wavevector that is identical for waves 
in different layers.— If we further assume that an oblique plane wave is incident from a 
dielectric medium with relative permitivity ei r , then k x = w^/ioeoeir sin(#), which can be 
substituted into Eq. §5§ to give, 

e{uj) — €i r sin 2 (#) 



dE x 
dy 



-iuno 7-^ H z . (6) 

Equations ([3]) and ([6]) constitute a system of equations for ID TM wave propagation across 
the interfaces between the media. To translate those equations into FDTD equations, Jiang 
et al. introduced a convenient method to overcome the difficulty of time-domain convolution 
between the term in the square bracket and H z in Eq. OH]). In this method, 29 Eq. (jH]) is first 
split into a pair of equations as, 

dE x 



dy 



-iuj^ H' z , (7) 
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H * ~ H » (8) 

Equations (J3J), (J7J) and jSJ) then lead to the following set of FDTD equations, 
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In obtaining Eq. (|T2|) . we have multiplied both sides of Eq. ([8]) by e(u) and Fourier 
transformed the result into the time domain. We have also made the assumption that a 
Drude model is used, e m = e(oo) — u 2 D /(u) 2 + iY D u). The updating coefficients in Eq. ( Ti"2"j) 

are 



'Jy3 



J y5 



'Jy4 



JyG 



7 [e r - 



€i r Sill 



2/ 



(13) 



in vacuum (e r = 1) and dielectric media (constant e r ), and 
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in Drude media. Here, we note the similarity between the updates of the (H', H) pair and 
the (P, D) pair in the UPML formulation, which results from the fact that both pairs involves 
updating an auxiliary variable before the treatment of the material dispersiveness. In the 
case that e(w) contains a linear sum of different types of poles (e.g., to accurately describe 
metals near inter-band transition energies^), direct Fourier transform may not be as efficient 
because of higher-order derivatives with respect to time. For a systematic treatment of this 
situation, interested readers are referred to Appendix O The updating coefficients in Eqs. 
((9]) to ( TTTT) are identical to those in Eqs. ( 1B1I) . ( 1B2I) and ( 1B6j) . which are given in Eqs. ( 1B81 
IB12j) . We note that the updating coefficients corresponding to Berenger's PML formulation 
can be used here provided that the two end media in the layers are vacuum. In the case 
of non-vacuum semi-infinite media at the two ends, ID UPML can be used to effectively 
absorb the outgoing waves, for example, equations similar to Eqs. (1B13HB151 IB 191 IB20I) can 
be used by setting n x = 1 and a x = in Eqs. ( 1B25I) and ( 1B27I) . 



B. Simulation setup and flow chart 



The main panel of Fig. [T] illustrates the geometry of the FDTD simulation region. The 
layered media are denoted by e ir , e 2r , etc. and are stacked along the y direction. The thick, 
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dashed (thin, dotted) lines denote the TF/SF boundaries along which the incident if- field 
(S-field) is calculated. Incident field alignments on the boundaries are shown more explicitly 
in the zoom-in panels to the left and below the main panel. In this work, we assume that 
the oblique incidence field is introduced from the lower left corner (crossing point between 
lines b and / in the main panel of Fig. [T]) with incident angle 9 to the normal of the media 
interfaces (y direction). We further assume that the two end media in the layers are vacuum. 
Consequently e\ r = 1, so that field propagation along the horizontal boundaries e through h 
can be readily calculated by a delay of the free-space propagation time. In addition, the ID 
field propagation along the vertical lines can be terminated by Berenger's PML formulation. 
The perfectly matched layers absorbing boundaries are not shown in Fig. [TJ They will 
be further illustrated and explained when we consider specific examples in Section IIHI The 
lower left panel in Fig. [1] shows the field alignment along line a for the ID wave propagation. 
The same setup applies to lines b, c and d. Importantly, a ID total field / scattered field 
approach is used here (the boundary points are highlighted in the shaded rectangle) because 
we must allow the wave from the multiple interface reflection to exit the ID simulation and 
be absorbed at the bottom on the ID simulation line.— 

Our simulation follows the flow chart shown in Fig. [2j The procedures belonging to ID 
and 2D field updates are highlighted in the shaded rounded rectangles. In each iteration, 
the code updates the ID i?-field, 2D .E-field, ID if-field, and 2D if-field in a sequence. The 
order of ID field storage and its matching to 2D simulation are important to ensure correct 
implementation of the 2D TF/SF scheme. Before updating the ID field, the code needs to 
store at each time instant the ID field values at the crossing points between line a and lines 
e, /, g, and h. 

The field matching at the TF/SF boundary is performed differently in accordance with 
the different updating schemes introduced in Section Hi Al In the ADE approach, the TF/SF 
boundary matching equations on lines e, /, b read, 

rpscalra+l rpscalri+l rrmclw+1/2 /i r\ 

^x \i+l/2,jl ~ x \i+l/2,jl a x2-n- zy |j + i/2,jl+l/2> \ LO ) 

rrtotin+3/2 _ rxtotin+3/2 pscam.+ l 

n zy li+l/2j'l+l/2 ~~ n zy lj+l/2j'l+l/2 "ifi^x H+1/2JV \ LU ) 
rrtotin+3/2 rrtotin+3/2 pscain+l 

n zx lil+l/2,j+l/2 n zx lil+l/2J+l/2 v x2^ y Uij+1/2' V 1 '/ 

These updates are performed immediately after Eqs. (IBip . (1B6I) . and (1B5I) . For the ^-field 



update on lines V and c', because J y depends on the updated value of E y in Eq. (IB4I) . the 



Ey boundary matching is performed in between Eqs. ( TB3I) and (1B4j) . for example, on line b', 

rpscaln+l pscalrt+l rrinclW+1/2 /io\ 

a y lilJ+1/2 ~ ^y lilJ+1/2 a y1 n zx lii+l/2,j+l/2- l ic V 

In the UPML formulation, the TF/SF boundary matching is carried out immediately 
after the P and B updates (before updating D and H) in Eqs. (1B13HB201) . for example, the 
updates on lines b and b' read 

pscaln+l _ psca|n+l rrinc|™+l/2 /iq\ 

y Ul,j+l/2 ~~ r y Ulj+1/2 a y^ n z lii+l/2,j+l/2' v iJ / 

ptot|«+3/2 _ R tot 1^+3/2 n pscaln+l (OfW 

D z l«l+l/2J+l/2 ~~ °z lil+l/2,i+l/2 P^y U1J+1/2- ^ ZU J 

Because the above updates are performed between the updates of P and D or B and H, 
they are indicated in the flow chart (Fig. [2]) by the upward arrows on the right. We note 
that if the same type of PML absorbing boundary condition is used to terminate both the 
ID and the 2D field propagation, one can allow them to have the same updating coefficients 
in the PML region and therefore remove the procedures of saving and matching the field 
components on lines g and h [-Ebot and Hbot}-— This particular setup is useful in the 
description of a very thick bottom layer (semi- infinite in the positive y direction). 

In the case of normal incidence, the code simplifies in two ways. First, in Eq. (JSJ), H = H', 
and therefore Eq. (JT2l) is removed from the ID if-field update procedure. Second, it is not 
necessary to store and interpolate the field values Ebot, -^top, Hbot, an d -f^rop, because 
field excitation is synchronized along lines e, /, g and h, respectively. 

The incident field values on lines a, b, c and d are calculated from Eqs. (l9 tiT2|) us- 
ing a ID TF/SF scheme that allows fields reflected from the interfaces to exit the ID 
simulation domain. Based on the geometry shown in the lower left panel in Fig. [1] 
we assume that the incoming if-field with time-dependence f(t) excites the ID field 
at point (il — 1/2, jl — 3/2). Paired with this excitation is an E-field of the form 



g(t) = — V/U /e /(t + Axcos(6 l )/2c) cos(6 ) ), exciting the ID field at point (il — 1/2, jl — 2). 
For example, on line a in Fig. [TJ the ID TF/SF boundary matching equations read, 

r-iscain+l r-iscam+l _ fin+1/2 /Q1 \ 

K \jl-2 = h a ljl-2 — a x-2j | , l/l) 

rr/totin+3/2 oTtot | n+3/2 j „\n+l / 'o<-A 

U a ljl-3/2 - U a ljl-3/2 _ °lfl9\ ■ \ ZZ ) 

The values of /|"+V2 an d 

g\ n are calculated from the known expressions of fit) and g(t) 



using time-domain interpolation when necessary. These values are stored at each instant 
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to generate the excitation fields for lines b, c, and d by introducing a time delay tdeiay = 
NAxsm(9)/c. The field values on the horizontal lines e, /, g and h are obtained in a 
similar fashion. One can also store the field values at each point on line a, and save the 
computation along lines b, c and d by introducing a proper time delay. This scheme reduces 
the computation time for the cost of larger memory requirement. Finally, the 2D TF/SF 
boundary values [E y nc ] along lines b' (d) are readily calculated from the if -field values on 
lines a and b (c and d) using Eqs. flB3j) and (IB4I) . In addition, we note that the excitation 
and PML absorbing boundary conditions are enforced on H' in the ID field updates. 

Several practical issues should be considered. First, in vacuum, the projection of the 
phase velocity of the oblique incident field on the y-axis is c/ cos(#). As the incident angle 
9 increases, the phase velocity can be very large and cause numerical instability if a fixed 
Courant criterium is enforced (e.g., At = Ax /2c). Based on this observation, we vary the 
Courant number S = cAt/Ax to ensure stability. When the incident angle is small, a small 
Courant number is used to ensure resolution of the time domain interpolation along the 
horizontal boundaries. In our simulation, the same Courant number is used for both ID and 
2D wave propagations, while an interpolation scheme to match different Courant numbers 
in ID and 2D wave propagations is explained in Ref. 28|. Second, as the dielectric function 
is discontinuous across the interface between layers of different media, we have used an 
average dielectric function for updating the fields at the interface.— 1 ^ For example, in the 
left panel of Fig. [Q we use the dielectric function e e g = (ei r + €2 r )/2 for the .E^-field updates 
at the interface. In Section IHIl we will show that this scheme leads to faster convergence 
and/or higher accuracy as compared to the standard step-like change of e. Finally, we use 
a Gaussian ramping in the hard source time-response in Eqs. (12 ip and (|22|) to slowly ramp 
the field to continuous wave so as to avoid high-frequency component leakage out of the 
TF/SF domain. Specifically, f(t) = exp(-(£ - T dealy ) 2 / t%) sin(wt), T dday = 30 fs, r = 10 fs, 
for the ramping phase < t < Tdeiay 



C. Numerical tests 



To test the accuracy of the TF/SF scheme, we compare our simulation results to analytical 
results by considering the oblique TM wave incident upon a slab sandwiched between two 
vacuum media. The analytical results for the reflection and transmission coefficients are 



9 



given by2 

ru + r 23 e 2 ' 7 
r ~ l + r 12 r 23 e 2 ^' 

1 + r 12 r 23 e 2 ^ ' 

For TM wave, 

cos(0j)/rij — cos(9j)/rij 

13 cos(9i) I rii + cos(0j) / % ' 

2cos(gj)/rai 

lJ cos(6 l i) /rij + cos(0j) / rij ' 
ft ' 

sin(6 l j) = — sin(^j), 

7 = —n 2 hcos(9 2 ). 

c 

Here, and t^- denote, respectively, the reflection and transmission coefficients at the 
interface between media i and j, and n« = ^e^w) denotes the refractive index of media 
i. We assume that medium 2 is a slab of thickness h. The waves at the input side of the 
slab, where the incident and reflected waves propagate, and at the output side, where the 
transmitted wave propagates, can then be expressed as, 

^input = exp ^ kxX + ky y)] + rexp[i(k x x - h y y)\, (29) 
^output = texp[i{k x x + k y y)}. (30) 

The above expressions indicate that the maximum field amplitude on the input and output 
sides are 1 + \r\ and \t\, respectively. These quantities can be obtained along a y-direction 
detection line in the TF/SF scheme for layered media without placing any scatterer inside the 
TF region. Using this scheme, we also test the leakage, defined as the ratio of the maximum 
field magnitude in the scattered field region to the maximum field magnitude in the total 
field region: leakage = max|^ sca, |"-|/max|^' tat |"-|, where ip refers to E x , or E y , or H z . In 
the ideal case, leakage = 0, whereas in practice, leakage < 10 -2 (or —40 dB) is desirable.— 
Specific numerical examples of the tests are provided in Section III II where "leakage" refers 
to the largest leakage among E x , or E y , or H z . In addition, we have tested the accuracy 
of the wave propagation in the layered media by inspecting the x and y projections of the 
wavelength [where, for instance, the x-projection is 2ir/k x and k x = ksm(6)}. In the case 
of a dispersive slab, we have also tested the skin depth (the distance where the field decays 
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(23) 
(24) 

(25) 

(26) 
(27) 
(28) 



to e _1 of its value at the surface, ca. 30 nm for the Drude model and parameters in our 
calculation), by considering a slab with thickness larger than 300 nm. These tests all show 
an error within 5% compared to analytical results. 

III. NUMERICAL EXAMPLES AND DISCUSSIONS 

To illustrate the generality of our formulation, we first consider the simple case of plane 
wave propagation in vacuum, illustrated in Fig. [3J In panels (a-c), the plane wave (wave- 
length 400 nm) is injected from the lower left corner into the TF region (bounded by the 
thick, dashed lines) with incident angle 6 = 65°. In panel (d), the plane wave propagates 
in the positive y direction. It is shown that as the field penetrates into the Berenger PML 
located at the top of the simulation domain, it is efficiently absorbed. Negligible leakage 
is introduced at the PML boundary as the ID field updating equations acquire the same 
coefficients as the 2D equations [see Section [IT] . The dashed oval in panel (a) indicates 
considerable leakage (3.292 x 10~ 2 ) outside the TF region because in this case the incident 
continuous wave (cw) field is turned on instantaneously. Consequently, the high frequency 
components in the leading wave front are not well matched at the TF/SF boundary, resulting 
in the leakage. As shown by panels (b) and (c), the leakage can be reduced by one order of 
magnitude by slow (Gaussian) ramping of the incident field to steady-state cw oscillations. 
In light of this, hereafter we use Gaussian ramping prior to cw in the excitation hard source 
f(t) and g{t). The maximum leakage in the calculations of panels (b-d) is 1.367 x 10~ 3 and 
the relative error in the vacuum wave impedance (Z = \J /UoAo) is 0.41%. 

As a second example, we study a plane wave obliquely incident on a dielectric slab.— In 
Fig. HI we plot snapshots of the magnetic field of a plane wave (wavelength 400 nm) incident 
at an angle 9 = 45° on a 900-nm thick dielectric slab (dielectric constant e r = 11.7). In both 
panels, the solid rectangle indicates the location of the slab, while the thick, dashed rectangle 
shows the TF/SF boundary. The plane wave is injected from the lower left corner and first 
impinges on the lower vacuum/dielectric interface. In panel (a), we observe the interference 
patterns of the reflected wave with the incident wave below the lower vacuum/dielectric 
interface while the refracted wave front propagates in the slab. The faint wave front in the 
dielectric slab is due to the slow Gaussian ramping of the incident field. After a steady state 
is established [Fig. SJb)], the magnetic field pattern clearly reveals the interference between 
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the reflected and incident waves, the interference within the dielectric slab, and the final 
transmission through the slab. In Fig. IU(b), it is observed that the final transmitted wave 
maintains the same propagation direction as the incident wave (45° to +y direction) because 
the media below and above the slab are both vacuum. 

We proceed to examine the convergence of the magnitude of the reflection (r) and trans- 
mission (t) coefficients to the analytical results given by Eqs. (|23l) and (1241) . In Fig. [5j we 
plot the relative errors in (a) |r| and (b) \t\ with respect to analytical results as a function 
of the mesh size Ax. The red, solid (blue, dashed) curve in Fig. [5] shows the convergence 
result without (with) the interface averaging of the dielectric constant. From the compar- 
ison, it is clear that calculations with interface correction lead to uniformly smaller error 
than that without the interface correction. The slope of each line in the log-log plot ob- 
tained by the least-square fit indicates that second order accuracy of Yee's algorithm is 
maintained with the interface correction, while the accuracy degrades to first order with- 
out the interface correction. Similar effects have been reported in previous studies on the 
accuracy of FDTD results with dielectric interfaces,— 1 ^ while here we observe such effects 
within the TF/SF formulation in the context of layered media. We note that the interface 
correction scheme does not entail additional computational and memory requirements and is 
thus always recommended. In Fig. [5](c), we show that the maximum leakage with interface 
averaging is uniformly smaller than that without the interface averaging for different mesh 
sizes. Throughout, the maximum leakage is below 2.0 x 1CT 3 , substantiating our confidence 
in the TF/SF scheme.— 

To further test the accuracy of the dielectric slab reflection and transmission upon oblique 
incidence of a plane wave, we compare the analytical results with FDTD calculated results 
at different incident wavelengths in Table |T] and at different incident angles in Table [TTJ As 
shown, the relative error (given in parentheses) is uniformly below 5%, except for incident 
angle 6 = 50°, where |r| is below 0.01. It is interesting to note that the relative error 
diminishes with increasing wavelength, while a non-monotonic trend is seen in the errors of 
both reflection and transmission coefficients for an increasing incident angle. 

Next we apply the TF/SF method to study the reflection and refraction of a plane wave 
obliquely incident upon a dispersive metal slab. Snapshots of the magnetic field are shown 
in Fig. [6] as the plane wave passes through the metal slab. Specifically, we consider an 
incident plane wave with wavelength 400 nm and 9 = 45° injected from the lower left 
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corner of the TF region upon an 80 nm thick dispersive metal slab described by the Drude 
model e m = e(oo) - u 2 D /(u 2 + iT D u), with e(oo) = 7.0246, u D = 1.5713 x 10 16 rad/s, and 
Tjj = 1.4003 x 10 14 rad/s. This set of parameters is optimized to fit the dielectric data 
reported in Ref. for bulk silver in the spectral range from 330 to 500 nm. Figures M (a) 
and (b) illustrate the magnetic field distribution before and after reaching a steady state, 
respectively. In Fig. [6] (b), the large curvatures at the interference minima between the 
incident and reflected fields below the lower interface indicate a large reflection coefficient 
(> 0.9). Inside the metal, because of the complex dielectric function of the slab, the wave 
front is no longer a plane wave, as is clearly discernable in Fig. El However, the final 
transmitted wave exiting from the upper interface recovers a plane wave front and the same 
propagation constant as the incident wave, because the media below and above the dispersive 
slab are both vacuum. From Figs. S(b) and [6(b), it is seen that the periodicity in the x 
direction of the fields below, inside, and above the slab is the same. By further observing 
the field propagation after reaching the steady state in both cases (not shown), it is clear 
that the phase of the total field in the x direction is matched. This observation confirms the 
phase matching condition parallel to the interface (same k x across the interfaces), which is 
critical to the derivation of the ID field propagation, Eq. ([6]). 

To examine the convergence of our results in the case of the metal slab, we use the 
same incident field condition as that in Fig. [6] and plot the relative error of the reflection 
and transmission coefficients as a function of the mesh size Ax in Figs. [7] (a) and (b), 
respectively. It is seen that the results with interface averaging (blue, dashed curves) of 
the dielectric function yield uniformly lower error than the results without the interface 
averaging (red, solid curves). A first-order power law is seen in the error of the transmission 
coefficient as a function of the mesh size without interface averaging, all other errors are 
near and below 10~ 3 , illustrating the convergence of the FDTD results. FDTD simulations 
on similar dispersive systems have been reported by Mohammodi et al., who suggested that 
the dispersive contour-path method is able to achieve smaller error even for a relatively 
large step-size (Ax).— Fig. [7(c) shows that the leakage decreases with a decreasing mesh 
size, albeit in this case the leakage with the interface averaging of the dielectric function is 
slightly larger than that without the averaging [cf. Fig. G3^c)]. 

In Tables [TTT1 and ITVt we compare between the analytical and the FDTD calculated reflec- 
tion and transmission coefficient magnitudes at various incident wavelengths and incident 
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angles for the metal slab studied in Fig. [6j The FDTD results are obtained after steady 
state is reached under cw incident plane wave illumination. In the frequency domain, this 
corresponds to a fixed incident wavelength, and the Drude model provides a constant com- 
plex value of dielectric function, which can be used in Eqs. f l2"3"j) and (|24l) to obtain the 
reflection and transmission coefficients. In Table 1111} we list the free space wavelength in 
the 350 to 500 nm range, to which the fitted Drude model is applicable. The small relative 
errors (< 2.5%) shown in the parentheses in Tables [TTT1 and HVl illustrates the reliability of 
our calculations using the TF/SF formulation in the case of layered dispersive media. The 
maximum leakage found in obtaining the data in Tables ITTT1 and HVl is 1.201 x 10~ 2 , which 
occurs at 6 = 80°. 

Panels in the left column of Fig. |H] illustrate snapshots of the magnetic field as the wave 
propagates through two-layer media consisting of a lower layer of 80-nm thick dispersive 
material and an upper layer of 100-nm thick dielectric material under oblique plane wave 
incidence. The material parameters are given in the caption of Fig. [8l In these panels, 
the solid horizontal lines define the boundaries between different layers, which are extended 
into the UPML in the x direction. The dashed box denotes the TF/SF boundary. The 
incident plane wave with wavelength 400 nm and 9 = 65° is injected from the lower left 
corner of the TF region. The magnetic field snapshots in the first column of Fig. [8] show 
that Snell's law is obeyed when the field passes through the two layers of materials. In 
particular, the propagation direction in the high-index dielectric material exhibits a smaller 
angle to the normal than the incident wave, whereas the final transmitted wave propagates 
along the direction of the incident wave. More importantly, we observe that the phase of 
the waves across the different layers is matched in the x direction after a steady state is 
established (bottom panel in the left column), which is again consistent with Eq. (Q. In 
this case, the magnitude of the reflection and transmission coefficients calculated by FDTD 
is |r| = 0.9524 and \t\ = 0.0866, respectively. The bottom panel in the left column also 
shows non-negligible leakage penetrating through the TF/SF box and propagating into the 
lower right corner of the simulation domain, nevertheless, the maximum value of the leakage 
in H z is 1.5461 x 10~ 4 , which is insignificant in practice. 

Panels in the second column of Fig. [8] are obtained under the same conditions as those in 
the first column except that a slit of 200 nm width (in the x direction) and 120 nm depth (in 
the y direction) is placed in the middle of the simulation domain. In the TF region, the slit 
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causes strong scattering of the injected plane wave, which results in the observed interference 
patterns. The slit introduces entirely new physics: outside the TF/SF boundary, the purely 
scattered wave distribution is reminiscent of a dipole radiation pattern. Closer inspection 
reveals that the field distribution is asymmetric with respect to the slit center (x = 600 
nm). The scattered field is strongest near the lower surface of the dispersive slab and to 
the right of the TF/SF box and weakest above the upper surface of the dielectric slab and 
to the left of the TF/SF box. The asymmetric angular distribution is a clear signature 
of the oblique incidence of the exciting plane wave. By enlarging the SF region size, we 
find that the purely scattered wave along the lower surface of the metal thin film consists 
mainly of surface plasmon polariton (SPP) waves propagating away from the slit. These 
are identified by their wavelength - 348 nm in the x direction compared with the analytical 
result for the wavelength of SPP at the interface between vacuum and metal, which is given 
by \ /Rey/e m /(l + e m ) = 348.25 nm. 

For the simulations in Fig. El we have updated the field at the horizontal and vertical 
interfaces using the averaging scheme discussed above, and have tested the convergence 
of the fields in the TF and SF regions with respect to mesh size (Ax), TF box size, and 
physical size of the simulation region. We note that UPML termination of the simulation 
domain is important because the scattered field due to the slit is significant. Our tests 
show that the maximum scattered field in the SF region is only one order of magnitude less 
than the maximum field in the TF region. Additionally, the UPML can effectively absorb 
the outgoing wave in the dispersive and dielectric layers. Furthermore, the boxed TF/SF 
boundary has advantage over the Il-shaped boundary considered previously, particularly 
when one is interested in the full angular distribution of the scattered field in the far-field 
zone. 

IV. CONCLUSIONS 

Using Maxwell's equations for the transverse magnetic wave, along with translational 
invariance and phase matching principles, we derived an equivalent one-dimensional wave 
propagation equation along the direction perpendicular to the interfaces between layered 
media. We then derived the corresponding finite-difference time-domain equations for lay- 
ered dielectric media and dispersive media with a Drude pole pair. To utilize these equations 
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for a plane wave with oblique incidence, we discussed the simulation setup and procedure 
in the framework of the total field / scattered field formulation with a special emphasis 
on techniques to match the fields at the total field / scattered field boundary. We have 
performed tests on vacuum propagation and on the reflection and refraction at a dielectric 
and a dispersive slab. Converged simulation results for various incident angles and wave- 
lengths reveal that the errors in the reflection and refraction coefficients are uniformly below 
5% compared to analytic results. The numerical example of scattering at a nano-scale (sub- 
wavelength) slit in a dispersive medium invites interesting applications of our formulation to 
time-dependent studies of electromagnetic wave scattering at surface or embedded scatters 
in dispersive media, for example, the coupling of incident oblique plane wave into surface 
plasmon polaritons. For this purpose, the developed method offers the flexibility of choos- 
ing the total field region inside which the near-field exhibits interference pattern between 
incident and scattered fields, while outside which the scattered far-field can be detected at 
all angles. 
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Appendix A: TE mode wave propagation 

Maxwell's equations in 2D in the frequency domain for the TE mode read, 



dH, 



y 



dH x 



(Al) 



dx 



dy 
dE z 



Z1 



dy 
dE z 



XI 



(A2) 



dx 



(A3) 
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Substituting Eq. ([A3]) into Eq. flAl]) yields, 

dH 

= iue [e(w) - e lr sin 2 (fl)] E z , (A4) 

where e\ r denotes the relative permitivity of the first medium (see Fig. [1]). Equations ( 1A2|) 
and (1A3I) are used for ID TE mode wave propagation. These equations can be readily 
solved using the same FDTD procedure as for Eqs. (j2]) and ([3]). The time-domain solution 
is facilitated by the fact that material dispersiveness introduces a factor [e(u) — t\ r sin 2 (6')] 
in Eq. (IA4[) . whereas in Eq. (E]) it introduces a factor [e(cj) — t\ T sin 2 (6*)] /e(w), which entails 
more difficulty for the FDTD solution.— The simulation setup and flow chart in Section [TT] 
can be used for the TE mode by exchanging the roles of the E and H fields. 



Appendix B: FDTD equations for 2D TM mode wave propagation 

The FDTD equations based on the auxiliary differential equation (ADE) approach read, 

P \n+l _ n rp in , / tt \n+l/2 jj in+1/2 n. j ,n+l/2 / R ..\ 

^H+l/2,3 ~ U xl nj x\i+l/2,j ^ u x2\ rl zy\i+l/2,j+l/2 n *V I i+l/2,j - 1/2 / ~^ a xZ<Jx \ i+ y 2 J \ DL J 

j \n+3/2 _ j in+1/2 . F m+1 m 9 \ 
t7a; lj+l/2j ~~ x ^ x H-{-l/2,j x ^ x H-\-l/2,j J y DZ ) 

n _L n (H r +1 / 2 - U \ n+1 / 2 

iJ+1/2 u l/2l- n 2x|j +1 /2 ) j + i/2 rl «*U-l/2, ; 



/? |n+l ej 
2/U,j+l/2 Uyl-^y 



j m+3/2 



1+3/2 - „ /■ r +1/2 4- „ p i n+1 

j+1/2 "1/4^^1^+1/2 "i" "'l/5- L 'j/lij+.i/2J 



-1/2 \ I t in+1/2 m „\ 

1/2J+1/2J «^3^y 1 1^+1/2 A- 00 ^ 



tt in+3/2 _ L rr in+1/2 . , / p m+1 p m+1 \ 

- n 2x|j+ 1 /2,j+l/2 ~~ " :cl - n ^U+l/2,j+l/2 "r ^ , 2l- C/ ?/U+l,j+l/2 • C/ l/liJ+l/2>'> 

in+3/2 _ , rr in+1/2 . , / p m+1 _ rp m+1 \ 

■y\i+l/2,j+l/2 u y 11Iz y\i+l/2,j+l/2 < u 2/2^:r|i+l/2,f+l ^U+l/2,j7' 



^2 = + H; 



zy 



(B4) 
(B5) 
(B6) 
(B7) 



The coefficients in the i?-field updating equations are medium dependent; specifically, in 
vacuum (e r = 1) and dielectric media (constant e r ),— 



Uxl — O-yl — 1; 

a X 2 = ~a y2 = At/(e e r Ax), 

dx3 — &x4 — a x5 — a y3 = a y4 



(B8) 
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in Drude media, e m = e(oo) — uj 2 d /{uj 2 + iTou) 
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Q'xl — a yl — 1) 

a X 2 = -a V 2 = At/ [e e(oo)Ax] , 
a X 3 = a V 3 = -At/ [e e(oo)] , 
a xi = a y4 = (1 - r D At)/(l + T D At), 
a x5 = a y s = e uj 2 D At/(l + T D At), 



(B9) 



and in the PML region, 



da* = exp(-o-^At/e ), 
a X 2 = [1 - exp(-cr 2/ At/e )] /(Axa y ), 

«a;3 — O x 4 = <Xj;5 = 0; 

a y i = exp(-a x At/e ), 
a y2 = ~ [1 - exp(-a x .At/e )] /(Axcr x ), 

a y3 — a yi — 0-y5 = 0. 

The coefficients in the if-field updating equations outside the PML regions are, 



(BIO) 



b x i 



b X 2 = -b y 

whereas in the PML regions they read, 



-At/(p Ax) 



(Bll) 



6a,! = exp(-o-*At//io), 

b x2 = - [1 - exp(-<7*At/// )] /(Axa*) 

b y i = exp(-cr*At/ / u ), 

b y2 =[l- exp(-^At//i Q )] /(Axa*). 



(B12) 



Here, we assume a polynomial grading of the PML parameters:— a X)V = eocr* //i = 
a m (p/S) m , where a m is the maximum conductance in the PML, p is the distance into the 
PML, and 5 is the thickness of the PML region. In this paper, we use a power m = 4 and 
5 = 20Ax. <j m is optimized to give a maximum reflection error on the order of 10 -7 . 
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The FDTD equations based on the UPML formulation read,— 



B- 



H- 
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The coefficients in the .E-fleld updating equations in all media are, 

a x \ = —dyi = At/ Ax. 
In vacuum (e r = 1) and dielectric media (constant e r ), 

■ ax2 = ail = at5 = aM = ai2 = atl = ayb = atS = , 



whereas in Drude media, 



0Lx2 


— a y% 


at X 3 


= oi y3 


®x4 




&x5 


= a yb 
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= OiyQ 



4e(oo) 

2e(oo)+e(oo)T D At+ui 2 D At 2 ' 
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2+r D At 
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2e(oo)+e(oo)r D At+w 2 D At 2 ' 

2-r D At 

2e(oo)+e(oo)r D At+w|,At 2 ' 



outside the UPML regions, 
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and in the UPML regions, 



Ot x 7 
®x9 



2eOKy—(TyAt 
2eOKy+(TyAt ' 

1 

2e^K y +eocr y At ' 

a x At + 2e n x , 



a x io = o- x At - 2e K x , 



a 



2eQK x —cr x At 
2eQK x +a x At ' 

1 

2ef ) K x +tocr x At ' 

a y At + 2e n 



(B25) 



«zio = c^At - 2e K y . 
The coefficients in the if-field updating equations outside the UPML region are, 



fa 


= & = 1, 


< ft 


= At/ Ax 


/4 


= 1Mb 



(B26) 



and in the UPML regions, 



fa 
fa 
fa 
fa 



2eQK x —<j x At 
2eoK x +a x At ' 

2e At 

(2e K x +a x At)Ax' 

2^0 — <7yAt 
2eo Ky-\-CTy At ' 

2€Q 



(B27) 



(2eoKj/+fvAt)/Lto ' 

Here, we assume a polynomial grading of the PML parameters,— a x>y 



{p/5) m and 



K x,y = 1 + (t m — where <r m and K m denote the maxima of the UPML parameters 

p is distance into the PML, and 5 is the thickness of the PML. In this paper, we use power 
m = 4, 5 = 20 Ax, K m = 1, and a m is optimized to give a maximum reflection error on the 
order of 10~ 7 for a simulation region consisting of vacuum and on the order of 10~ 4 for a 
simulation region consisting of Drude dispersive media. 
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Appendix C: Systematic solution of one-dimensional wave propagation in the TM 
mode 



In this appendix we provide a systematic solution for Eq. (JH]) when < sin(6 l ) < 1 and 
e(u) consists of a linear superposition of Debye, Drude, and Lorentz types of poles. In this 
case, we first rearrange Eq. (JS} as 



e lr sm 2 (8)H z = e(Lo)(H z -H' z ) } 



where 



e(u) = e(oo) + 



(CI) 



(C2) 



an 



d2a 



6i(u)) = < 



l-iu)T DB ' 
^DR 



d 2 -uj'j +i2F L u> ' 



for a Debye pole 
for Drude pole pairs 
for Lorentz pole pairs. 



(C3) 



We introduce auxiliary variables Ki to rewrite Eq. ( 1C1I) as a system of equations 



e lr sm 2 (9)H z 

K 



% 

e(oo)(H z -H> z ). 
e r (uj)(H z -H>). 



(C4) 

(C5) 
(C6) 



Equations ( 1C4I) and ( 1C5I) correspond to the set of FDTD equations, 

e lr sin 2 (^)^r +3 / 2 - K \ n+3 / 2 - J2^\ n+3/2 = 0, 

i 

e(oo)H z \ n+3 / 2 - K r +3 / 2 = e(oo)^| n+3 / 2 . 



(C7) 
(C8) 



The translation of Eq. ( 1C6I) into a set of FDTD equations depends on the type of pole(s) 
considered [see Eq. ( 1C3I) ]. For a single Debye pole (Ki = Kbb), 

At DB AtH z \ n+3 ' 2 - (2t db + At)K DB \ n+3 ' 2 
= (At-2r DB )K DB \ n+1 / 2 

+Ae DB At(H' z \ n+3/2 + H' z \ n+1/2 ~ H z \ n+l ' 2 ). (C9) 
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For a Drude pole pair (Ki = Kdr), 

u 2 DR At 2 H z \^ 2 - (T DR At + 2)K DR \^ 2 
= -AK DR \ n+l ^ 2 + (2 - V DR At)K DR \ n - x ' 2 

+u 2 DR At 2 (H' z \ n+3/2 + H' z \ n ~ 1/2 - H^- 1 ' 2 ). (CIO) 
For a Lorentz pole pair (K^ = K^), 

Ae L co 2 L At 2 H z \ n+3/2 - (cu 2 L At 2 + 2T L At + 2)K L \ n+3/2 
= -AK L \ n+1/2 + (u 2 L At 2 - 2r L At + 2)K L \ n - 1 ' 2 

+Ae L u 2 L At 2 (H' z \ n+V2 + H' z \ n ~ 1/2 - H z \ n ~ l/2 ). (Cll) 

Equations (1C7P through (101 lj) form a linear system of equations for the unknowns H z \ n+3 ^ 2 
and Ki\ n+ ^ 2 (i = 0,1,...), which can be solved by existing numerical solvers for linear 
systems of equations. The solution is then used to replace Eq. (|T2l to proceed the ID 
wave propagation for the TM mode Comparing to a direct Fourier transform of Eq. (jBj), 
the above procedure only requires the storage of the quantities at the previous two time 
instances and thus avoids the complexity of numerical high-order derivatives with respect 
to time. This procedure can be extended systematically to multi-poles in the material 
dispersiveness, although it involves solving a linear system of equations. 
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FIGURES 



b' 




a 



FIG. 1: Simulation geometry: the layered media are distinguished by different shades and denoted 
by €i r , £2r, etc. Thick, dashed (thin, dotted) lines denote the boundaries to which the H— (E-) 
field is assigned. The left, lower left, and lower right panels show the specific field point assignment 
at the interface, along line a, and at the horizontal boundaries, respectively. The x-coordinates 
of lines a, b', 6, c, c', d are i\ — 1/2, i\, i\ + 1/2, 12 — 1/2, and ii + 1/2, respectively. The 
y-coordinates of lines e, /, g, and h are ji, j\ + 1/2, ji — 1/2, and j'2, respectively. 
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2D Wave Propagation 



FIG. 2: Simulation flow chart. Note that "TOP" ("BOT") refers to a lower (higher) y coordinate. 
For ID field updates an Auxiliary Differential Equation (ADE) approach is used, while 2D field 
updates are performed by either the ADE approach or the equations consistent with the Uniaxial 
Perfectly Matched Layers (UPML) formulation. 
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FIG. 3: (color online) In panels (a-c), the incident plane wave enters the TF region from the lower 
left corner of the TF/SF boundary with incident angle 9 = 65°. (a) Magnetic (H) field snapshot at 
3.00 fs for a plane wave propagating in vacuum. The dashed oval indicates the leakage outside the 
TF/SF boundary as a result of the instantaneous turn-on of the field. Panels (b-d) show .ff-field 
snapshots for a plane wave propagation with initial Gaussian ramping, (b) ff-field snapshot during 
ramping (at 7.34 fs) and (c) after steady state is established (at 66.71 fs). (d) -ff-field snapshot at 
90.06 fs for a plane wave propagating in the positive y direction. For all calculations, the incident 
wavelength is A = 400 nm (period T = 1.33 fs), and the steady-state amplitude of the incident 
magnetic field is 1 A/m. The mesh size is Ax = 2.5 nm, and the Courant number is S = 0.4. A 
log color scale (log 10 \H Z \) is used in all plots. The thick, dashed (thin, dotted) rectangle indicates 
the TF/SF (inner PML) boundary. 
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FIG. 4: (color online) Magnetic field snapshots of a plane wave obliquely incident on a dielectric 
slab (indicated by a solid rectangle). Reflection and refraction (a) at the lower interface (at 10.01 
fs) and (b) after steady state is established (at 100.07 fs). The incident plane wave enters the TF 
region from the lower left corner of the TF/SF boundary with incident angle 6 = 45°. The incident 
wavelength in vacuum is A = 400 nm (period T = 1.33 fs), and the steady-state amplitude of the 
incident magnetic field is 1 A/m in all calculations. The dielectric constant and thickness of the 
slab are e r = 11.7 and 900 nm, respectively. The media above and below the slab are vacuum, 
the mesh size is Ax = 1.0 nm, and the Courant number is S = 0.3. A log color scale is used in 
all plots. The thick, dashed (thin, dotted) rectangle indicates the TF/SF (inner PML) boundary. 
The slab does not penetrate into the PML region. 
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FIG. 5: (color online) Relative error in the magnitude of (a) the reflection coefficient r and (b) the 
transmission coefficient i as a function of the mesh size Ax. (c) Maximum leakage as a function of 
mesh size Ax. In all calculations, the Courant number is S = 0.3. In all figures, the red, solid (blue, 
dashed) curve shows the result without (with) the interface averaging of the dielectric constants. 
The parameters of the incident wave and the dielectric slab are as in the calculation leading to Fig. 

SI 
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FIG. 6: (color online) Magnetic field snapshots of a plane wave obliquely incident on a Drude metal 
slab (indicated by a solid rectangle). Reflection, refraction and transmission (a) before (at 2.50 
fs) and (b) after (at 60.04 fs) steady state is established. In all calculations, the incident plane 
wave enters the TF region from the lower left corner of the TF/SF boundary (dashed lines) with 
incident angle 9 = 45°. The incident wavelength in vacuum is A = 400 nm (period T = 1.33 fs), 
and the steady-state amplitude of the incident magnetic field is 1 A/m. The metal slab is 80 nm 
thick with Drude parameters: e(oo) = 7.0246, lu d = 1.5713 x 10 16 rad/s, and T D = 1.4003 x 10 14 
rad/s. The media above and below the slab are vacuum, the mesh size is Ax = 1.0 nm, and the 
Courant number is S = 0.3. A log color scale is used in all plots. The thick, dashed (thin, dotted) 
rectangle indicates the TF/SF (inner PML) boundary. The slab does not penetrate into the PML 
region. 
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FIG. 7: (color online) Relative error in the magnitude of the (a) the reflection coefficient r and 
(b) the transmission coefficient i as a function of the mesh size Ax. (c) Maximum leakage as 
a function of mesh size Ax. The Courant number is S = 0.3 in all calculations. In all figures, 
the red, solid (blue, dashed) curve shows the result without (with) the interface averaging of the 
dielectric constants. The parameters of the incident wave and the metal slab are the same as in 
the calculations leading to of Fig. [6l 
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FIG. 8: (color online) Magnetic field snapshots for a plane wave obliquely incident on two layers of 
materials. The interfaces between the layers and vacuum are indicated by solid horizontal lines. The 
lower layer is an 80 nm thick Drude metal with parameters: e(oo) = 7.0246, cod = 1-5713 x 10 16 
rad/s, and = 1.4003 x 10 14 rad/s. The upper layer is a 100 nm dielectric with dielectric 
constant e r = 11.7. The media below and above the two layers are vacuum. The panels in the 
right column illustrate the scattering due to a slit of width 200 nm and depth 120 nm in the same 
layered structure as in the left column. Rows 1, 2, and 3 show snapshots at 1.60, 3.20, 4.80 fs 
(before a steady state is established); Row 4 shows the snapshot at 60.04 fs (after a steady state 
is established). For all calculations, the incident plane wave enters the TF region from the lower 
left corner of the TF/SF boundary with with incident angle 9 = 65°. The incidence wavelength 
is A = 400 nm (period T = 1.33 fs), and the steady-state amplitude of the incident magnetic field 
is 1 A/m. The mesh size is Ax = 2.0 nm, and the Courant number is S = 0.3. A log color scale 
is used in all plots. The thick, dashed (thin, dotted) rectangle indicates the TF/SF (inner PML) 
boundary. The slabs are extended into the UPML region. 
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TABLES 



TABLE I: Comparison of the magnitude of the reflection (r) and transmission (t) coefficients 
between the analytical and numerical results for different incidence wavelengths (A). Superscript 
a denotes the analytical, and superscript n denotes the numerical results. The percentages in 
brackets denote the relative errors in the numerical results. The mesh size is Ax = 1 nm and the 
Courant number is S = 0.3. 



A (nm) \r a \ \r n \ (error) \t a \ \t n \ (error) 

300 0.2485 0.2600(4.63%) 0.9686 0.9658(0.29%) 

400 0.1898 0.1949(2.69%) 0.9818 0.9809(0.09%) 

500 0.1532 0.1558(1.70%) 0.9882 0.9878(0.04%) 

600 0.1282 0.1298(1.25%) 0.9912 0.9916(0.04%) 

700 0.6990 0.6988(0.03%) 0.7152 0.7154(0.03%) 

800 0.7172 0.7171(0.01%) 0.6969 0.6969(< 0.02%) 
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TABLE II: Comparison of the magnitude of the reflection (r) and transmission (t) coefficients 
between the analytical and numerical results for different incidence angles (8) and Courant numbers 
(S). Results are obtained for a plane wave with 400 nm wavelength incident upon a 900 nm thick 
dielectric slab (e r = 11.7). Superscript a denotes the analytical, and superscript n denotes the 
numerical results. The percentages in brackets denote the relative errors in the numerical results. 
The mesh size is Ax = 1 nm. 



6 (degree) 
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0.3 


0.8278 


0.8281(0.04%) 


0.5610 


0.5539(1.27%) 


10 


0.1 


0.8169 


0.8174(0.06%) 


0.5767 


0.5760(0.12%) 


20 


0.2 


0.7737 


0.7747(0.13%) 


0.6335 


0.6323(0.21%) 


30 


0.3 


0.6565 


0.6588(0.35%) 


0.7543 


0.7526(0.23%) 


40 


0.3 


0.3839 


0.3884(1.17%) 


0.9234 


0.9216(0.19%) 
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0.3 


0.0039 


0.0089(128%) 
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1.0000(< 0.01%) 


60 


0.3 
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0.1952(1.46%) 
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0.9808(0.06%) 
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0.2 


0.1152 


0.1142(0.87%) 


0.9933 


0.9935(0.02%) 


80 


0.1 


0.3395 


0.3367(0.82%) 


0.9406 


0.9470(0.68%) 



TABLE III: As in Table [I] for oblique incidence upon an 80 nm thick silver slab. The Dielectric 
function of silver is described by the Drude model, e m = e(oo) — oo^,/(oj 2 + T^w), with e(oo) = 
7.0246, lo d = 1.5713 x 10 16 rad/s, and T D = 1.4003 x 10 14 rad/s. The mesh size is Ax = 5 nm 
and the Courant number is S = 0.3, except for the A = 350 nm case, where Ax = 2 nm and the 
Courant number is S = 0.1. 



A (nm) \r a \ \r n \ (error) \t a \ \t n \ (error) 

350 0.8908 0.8908(< 0.01%) 0.2295 0.2295(< 0.04%) 

400 0.9503 0.9500(0.03%) 0.1231 0.1234(0.24%) 

450 0.9672 0.9672(< 0.01%) 0.0760 0.0761(0.13%) 

500 0.9743 0.9742(0.01%) 0.0532 0.0533(0.19%) 
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TABLE IV: As in Table [U for oblique incidence upon an 80 nm thick silver slab. The dielectric 
function of silver is described by the Drude model, e m = e(oo) — uJ^jiJ 1 + iTduj), with e(oo) = 
7.0246, uj d = 1.5713 x 10 16 rad/s, and T D = 1.4003 x 10 14 rad/s. The mesh size is Ax = 5 nm. 
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